home *** CD-ROM | disk | FTP | other *** search
/ The PC-SIG Library 10 / The PC-Sig Library - Shareware for the IBM PC and Compatibles (PC-SIG)(Tenth Edition Disks 1-2804)(1991).iso / PC_SIGCD / 09 / 2 / DISK0921.ZIP / KEPLER.BAS < prev    next >
BASIC Source File  |  1987-09-30  |  5KB  |  211 lines
  1. ' PROGRAM "KEPLER"
  2.  
  3. ' NUMERICAL SOLUTION OF KEPLER'S EQUATION FOR
  4. ' ELLIPTIC, PARABOLIC AND HYPERBOLIC ORBITS
  5.  
  6. ' PUBLIC DOMAIN
  7.  
  8. ' IBM-PC  << QUICKBASIC COMPILER VERSION 3.0 >>
  9.  
  10. DEFDBL A-Z
  11.  
  12. CONST PI=3.14159265D0
  13. CONST RTD=180/PI
  14. CONST DTR=PI/180
  15.  
  16. COMMON SHARED ECC,MANOM.RAD,EANOM.RAD
  17.  
  18. '**********************************************************
  19.  
  20. CLS
  21. PRINT
  22. PRINT
  23. PRINT "Program KEPLER"
  24. PRINT
  25. PRINT "Microsoft QuickBASIC Compiler"
  26. PRINT "(C) Copyright Microsoft Corp. 1982-1987"
  27. PRINT
  28. CALL KEYCHECK
  29.  
  30. CLS
  31. PRINT
  32. PRINT
  33. PRINT "Program introduction ( y = yes, n = no )"
  34. INPUT A$
  35. IF INSTR("yY",A$) THEN CALL INTRODUCTION
  36.  
  37. DO
  38.    CLS
  39.    PRINT
  40.    PRINT
  41.    PRINT TAB(15); "Kepler's Equation"
  42.    PRINT
  43.    PRINT
  44.    PRINT "Mean anomaly ( degrees [ 0 - 360 ] )"
  45.    INPUT MANOM.DEG
  46.    MANOM.RAD=DTR*MANOM.DEG
  47.    PRINT
  48.    PRINT "Orbital eccentricity ( non-dimensional )"
  49.    INPUT ECC
  50.  
  51.    IF ECC<>0D0 THEN
  52.       CALL KEPLER
  53.    ELSE
  54.       EANOM.RAD=MANOM.RAD
  55.    END IF
  56.  
  57.    ' PRINT RESULTS
  58.  
  59.    FORMAT$="####.######"
  60.  
  61.    CLS
  62.    PRINT
  63.    PRINT
  64.    PRINT TAB(22);"Kepler's Equation"
  65.    PRINT
  66.    PRINT
  67.    PRINT TAB(5);"Eccentricity       ( non-dimensional )";TAB(45);
  68.    PRINT USING FORMAT$;ECC
  69.    PRINT
  70.    PRINT TAB(5);"Mean anomaly       ( degrees )";TAB(45);
  71.    PRINT USING FORMAT$;MANOM.DEG
  72.    PRINT
  73.    IF ECC<1 THEN PRINT TAB(5);"Eccentric anomaly  ( degrees )";TAB(45);
  74.    IF ECC=1 THEN PRINT TAB(5);"Parabolic anomaly  ( degrees )";TAB(45);
  75.    IF ECC>1 THEN PRINT TAB(5);"Hyperbolic anomaly ( degrees )";TAB(45);
  76.    PRINT USING FORMAT$;EANOM.RAD*RTD
  77.    CALL KEYCHECK
  78.  
  79.    CLS
  80.    PRINT
  81.    PRINT
  82.    PRINT "Another selection ( y = yes, n = no )"
  83.    INPUT SELECTION$
  84. LOOP UNTIL INSTR("nN",SELECTION$)
  85.  
  86. END
  87.  
  88. '**********************************************************
  89.  
  90. SUB KEPLER STATIC
  91.  
  92.     ' KEPLER'S EQUATION SUBROUTINE
  93.  
  94.     EANOM.RAD=LOG(2*MANOM.RAD/ECC+1.8D0)
  95.     IF ECC<1D0 THEN CALL ESTIMATE
  96.  
  97.     DO
  98.        IF ECC<1D0 THEN L=ECC*SIN(EANOM.RAD)
  99.        IF ECC>=1D0 THEN L=.5D0*ECC*(EXP(EANOM.RAD)-EXP(-EANOM.RAD))
  100.  
  101.        IF ECC<1D0 THEN M=EANOM.RAD-L-MANOM.RAD
  102.        IF ECC>=1D0 THEN M=L-EANOM.RAD-MANOM.RAD
  103.  
  104.        IF ABS(M)<=1D-8 THEN EXIT DO
  105.  
  106.        IF ECC<1D0 THEN P=ECC*COS(EANOM.RAD)
  107.        IF ECC>=1D0 THEN P=.5D0*ECC*(EXP(EANOM.RAD)+EXP(-EANOM.RAD))
  108.        IF ECC<1D0 THEN Q=1D0-P
  109.        IF ECC>=1D0 THEN Q=P-1D0
  110.  
  111.        R=-M/Q
  112.        R=-M/(Q+.5D0*R*L)
  113.        R=-M/(Q+.5D0*R*L+R*R*P/6D0)
  114.        R=-M/(Q+.5D0*R*L+R*R*P/6D0-R*R*R*L/24D0)
  115.        EANOM.RAD=EANOM.RAD+R
  116.     LOOP
  117.  
  118. END SUB
  119.  
  120. '**********************************************************
  121.  
  122. SUB ESTIMATE STATIC
  123.  
  124.     ' EMPIRICAL ESTIMATE SUBROUTINE
  125.  
  126.     A=ECC^2
  127.     B=MANOM.RAD^2
  128.  
  129.     EANOM.RAD=A*(-.0579437D0*B+.02665324D0*MANOM.RAD+.05868658D0)
  130.     EANOM.RAD=EANOM.RAD+ECC*(-.19174923D0*B+.46905672D0*MANOM.RAD+.600725329D0)
  131.     EANOM.RAD=EANOM.RAD-.04987627D0*B+1.17552263D0*MANOM.RAD-.12324871D0
  132.  
  133. END SUB
  134.  
  135. '**********************************************************
  136.  
  137. SUB INTRODUCTION STATIC
  138.  
  139.     ' INTRODUCTION SUBROUTINE
  140.  
  141.     CLS
  142.     PRINT
  143.     PRINT TAB(10);"PROGRAM 'KEPLER' IS AN INTERACTIVE BASIC"
  144.     PRINT TAB(10);"PROGRAM WHICH CAN BE USED TO SOLVE"
  145.     PRINT TAB(10);"KEPLER'S EQUATION. THE ORBITAL MOTION"
  146.     PRINT TAB(10);"OF SATELLITES, COMETS AND ALL CELESTIAL"
  147.     PRINT TAB(10);"BODIES CAN BE DETERMINED FROM KEPLER'S"
  148.     PRINT TAB(10);"EQUATION. THIS EQUATION IS WRITTEN AS"
  149.     PRINT
  150.     PRINT TAB(10);"       M = E - e SIN E "
  151.     PRINT
  152.     PRINT TAB(10);"FOR ELLIPTIC ORBITS OR"
  153.     PRINT
  154.     PRINT TAB(10);"       M = e SINH F - F"
  155.     PRINT
  156.     PRINT TAB(10);"FOR PARABOLIC OR HYPERBOLIC ORBITS."
  157.     PRINT
  158.     PRINT TAB(10);"KEPLER'S EQUATION IS A TRANSCENDENTAL"
  159.     PRINT TAB(10);"EQUATION WHICH IS USUALLY SOLVED WITH"
  160.     PRINT TAB(10);"AN ITERATIVE NUMERICAL METHOD."
  161.     CALL KEYCHECK
  162.  
  163.     CLS
  164.     PRINT
  165.     PRINT TAB(10);"   IN BOTH OF THESE EQUATIONS, 'M'"
  166.     PRINT TAB(10);"IS CALLED THE 'MEAN ANOMALY' AND 'E'"
  167.     PRINT TAB(10);"IS CALLED THE 'ORBITAL ECCENTRICITY'."
  168.     PRINT TAB(10);"IN THE FIRST EQUATION, 'E' IS CALLED"
  169.     PRINT TAB(10);"THE 'ECCENTRIC ANOMALY' WHILE IN THE"
  170.     PRINT TAB(10);"SECOND EQUATION 'F' IS CALLED THE"
  171.     PRINT TAB(10);"'PARABOLIC' OR 'HYPERBOLIC' ANOMALY."
  172.     PRINT
  173.     PRINT TAB(10);"THE TRIG FUNCTION 'SINH' IN THE SECOND"
  174.     PRINT TAB(10);"EQUATION IS THE HYPERBOLIC SINE."
  175.     PRINT
  176.     PRINT TAB(10);"   THE TYPE OF ORBIT IS DETERMINED"
  177.     PRINT TAB(10);"BY THE ECCENTRICITY."
  178.     PRINT
  179.     PRINT TAB(10);"       0 < e < 1   ELLIPTIC"
  180.     PRINT TAB(10);"           e = 1   PARABOLIC"
  181.     PRINT TAB(10);"           e > 1   HYPERBOLIC"
  182.     CALL KEYCHECK
  183.  
  184.     CLS
  185.     PRINT
  186.     PRINT TAB(10);"THE ACTUAL POSITION OF A BODY IN ITS"
  187.     PRINT TAB(10);"ORBIT IS CALLED 'TRUE ANOMALY' AND"
  188.     PRINT TAB(10);"CAN BE DETERMINED FROM ECCENTRICITY"
  189.     PRINT TAB(10);"AND THE ECCENTRIC OR PARABOLIC OR"
  190.     PRINT TAB(10);"HYPERBOLIC ANOMALY."
  191.     CALL KEYCHECK
  192.  
  193. END SUB
  194.  
  195. '**********************************************************
  196.  
  197. SUB KEYCHECK STATIC
  198.  
  199.     ' CHECK USER RESPONSE SUBROUTINE
  200.  
  201.     PRINT
  202.     PRINT
  203.     PRINT TAB(15);"< press any key to continue >"
  204.  
  205.     A$=""
  206.     WHILE A$=""
  207.       A$=INKEY$
  208.     WEND
  209.  
  210. END SUB
  211.